** estimate effects by user age

***Inputs: 
* $Data/GameUserWPdaily_norm.dta
* $Data/maxmindirectionsmonth_zip260.dta

***Ouptus: 
* $Results/byage.ster
* $Results/agedetail_int.png
* $Results/agedetail_intabove5plays.png
* $Results/agedetail_intgames.png
* $Results/agescore.tex
* $Results/agescore_obs.png
* $Results/ageobs.tex
	
	use $Data/GameUserWPdaily_norm.dta, clear
	
	gen month = month(date)

		
	gen zip2 = int( zip/10)
	
	merge m:1 zip2 month using $Data/maxmindirectionsmonth_zip260.dta
	
	keep if _merge ==3
	
	
	gen monthyear = mofd(date)

	gen highpm =  weighted_dailypm > 25 & weighted_dailypm<. 

	
	
	egen cut_winddirection = cut(median_closestwind), at(0(60)360)
	replace cut_winddirection = 0 if median_closestwind ==360

	
	
	gen maxzipwind = (cut_winddirection ==maxdirection )
	replace maxzipwind = . if maxdirection==.
	
	keep if maxzipwind<. 
	
	
		egen cut_temp = cut(weighted_temperature), at(-40(10)40)

	
label variable highpm "Daily PM2.5 $>$ 25" 
	
	
 egen cut_age = cut(age_sig), at(20(10)81)

 replace cut_age = 70 if age> =80
 replace cut_age = 20 if age<20
	
	
	
	**** Results by age 
	
	ivreghdfe norm_score (c.highpm#i.cut_age    =c.maxzipwind#i.cut_age  ) weighted_wind_speed weighted_precipitation , absorb(anon_id monthyear  cut_temp) cluster(anon_id monthyear)


	eststo agedetail_int

	estimates save $Results/byage, replace
	
	
	
	
coefplot  (agedetail_int, keep(20.cut_age#c.highpm) ciopts(recast(rcap) lcolor(red%30))  mcolor(red) msymbol(D) msize(large) omitted baselevels ) ///
	( agedetail_int, keep(30.cut_age#c.highpm) ciopts(recast(rcap) lcolor(orange%30))  mcolor(orange) msymbol(D) msize(large) ) ///
	( agedetail_int , keep(40.cut_age#c.highpm) ciopts(recast(rcap) lcolor(gold%30)) mcolor(gold) msymbol(D) msize(large)  )  ///
	(agedetail_int,  keep(50.cut_age#c.highpm) ciopts(recast(rcap) lcolor(green%30)) mcolor(green) msymbol(D) msize(large) ) ///
	( agedetail_int, keep(60.cut_age#c.highpm) ciopts(recast(rcap) lcolor(blue%30)) mcolor(blue) msymbol(D) msize(large) ) ///
	( agedetail_int, keep(70.cut_age#c.highpm) ciopts(recast(rcap) lcolor("111 000 255 %30")) mcolor("111 000 255") msymbol(D) msize(large) ) ///	
	,  /// 
	drop(weighted*)  ///
			     vertical levels(95) yline(0)   legend(off)      	 ///
				 coeflabels( ///
                20.cut_age#c.highpm = "<30" ///
			  30.cut_age#c.highpm = "30" ///
			  40.cut_age#c.highpm = "40" ///
			  50.cut_age#c.highpm = "50" ///
			  60.cut_age#c.highpm = "60"  ///
			  70.cut_age#c.highpm = "{&ge}70" ///
		) ///
			 nooffsets 
			 
	graph export $Results/agedetail_int.png, replace
	
	
			
	
	
	**** Results by age above 5 play effects

	
	
	gen above5play = nth_play >= 5 & nth_play<.

	
	ivreghdfe norm_score (c.highpm#i.cut_age    =c.maxzipwind#i.cut_age  ) weighted_wind_speed weighted_precipitation , absorb(anon_id monthyear  cut_temp above5play) cluster(anon_id monthyear)


	eststo agedetail_int

coefplot  (agedetail_int, keep(20.cut_age#c.highpm) ciopts(recast(rcap) lcolor(red%30))  mcolor(red) msymbol(D) msize(large) omitted baselevels ) ///
	( agedetail_int, keep(30.cut_age#c.highpm) ciopts(recast(rcap) lcolor(orange%30))  mcolor(orange) msymbol(D) msize(large) ) ///
	( agedetail_int , keep(40.cut_age#c.highpm) ciopts(recast(rcap) lcolor(gold%30)) mcolor(gold) msymbol(D) msize(large)  )  ///
	(agedetail_int,  keep(50.cut_age#c.highpm) ciopts(recast(rcap) lcolor(green%30)) mcolor(green) msymbol(D) msize(large) ) ///
	( agedetail_int, keep(60.cut_age#c.highpm) ciopts(recast(rcap) lcolor(blue%30)) mcolor(blue) msymbol(D) msize(large) ) ///
	( agedetail_int, keep(70.cut_age#c.highpm) ciopts(recast(rcap) lcolor("111 000 255 %30")) mcolor("111 000 255") msymbol(D) msize(large) ) ///	
	,  /// 
	drop(weighted*)  ///
			     vertical levels(95) yline(0)   legend(off)      	 ///
				 coeflabels( ///
                20.cut_age#c.highpm = "<30" ///
			  30.cut_age#c.highpm = "30" ///
			  40.cut_age#c.highpm = "40" ///
			  50.cut_age#c.highpm = "50" ///
			  60.cut_age#c.highpm = "60"  ///
			  70.cut_age#c.highpm = "{&ge}70" ///
		) ///
			 nooffsets 
			 
	graph export $Results/agedetail_intabove5plays.png, replace
	
	
	
	
	
	
	**** Results by age with domain effects

	
	ivreghdfe norm_score (c.highpm#i.cut_age    =c.maxzipwind#i.cut_age  ) weighted_wind_speed weighted_precipitation , absorb(anon_id monthyear  cut_temp gamenum) cluster(anon_id monthyear)


	eststo agedetail_int

coefplot  (agedetail_int, keep(20.cut_age#c.highpm) ciopts(recast(rcap) lcolor(red%30))  mcolor(red) msymbol(D) msize(large) omitted baselevels ) ///
	( agedetail_int, keep(30.cut_age#c.highpm) ciopts(recast(rcap) lcolor(orange%30))  mcolor(orange) msymbol(D) msize(large) ) ///
	( agedetail_int , keep(40.cut_age#c.highpm) ciopts(recast(rcap) lcolor(gold%30)) mcolor(gold) msymbol(D) msize(large)  )  ///
	(agedetail_int,  keep(50.cut_age#c.highpm) ciopts(recast(rcap) lcolor(green%30)) mcolor(green) msymbol(D) msize(large) ) ///
	( agedetail_int, keep(60.cut_age#c.highpm) ciopts(recast(rcap) lcolor(blue%30)) mcolor(blue) msymbol(D) msize(large) ) ///
	( agedetail_int, keep(70.cut_age#c.highpm) ciopts(recast(rcap) lcolor("111 000 255 %30")) mcolor("111 000 255") msymbol(D) msize(large) ) ///	
	,  /// 
	drop(weighted*)  ///
			     vertical levels(95) yline(0)   legend(off)      	 ///
				 coeflabels( ///
                20.cut_age#c.highpm = "<30" ///
			  30.cut_age#c.highpm = "30" ///
			  40.cut_age#c.highpm = "40" ///
			  50.cut_age#c.highpm = "50" ///
			  60.cut_age#c.highpm = "60"  ///
			  70.cut_age#c.highpm = "{&ge}70" ///
		) ///
			 nooffsets 
			 
	graph export $Results/agedetail_intgames.png, replace
	
	
	
	
	
	
	
	
	
	
	
	
	
	******** Distribution of age 
	
	
	
		gen ones = 1
		
		
		
		gen agebins = 1 if cut_age==20
		replace agebins = 2 if cut_age==30
		replace agebins = 3 if cut_age==40
		replace agebins = 4 if cut_age==50
		replace agebins = 5 if cut_age==60
		replace agebins = 6 if cut_age==70
		
		
		label define bins 1 "$<$30" 2 "30"	 3 "40"	 4 "50"	 5 "60"	 6 "$\geq$>70"	

				label values agebins bins 

		
		reg norm_score i.agebins if maxzipwind<. 
		
		eststo agescore
		
		
		
	esttab agescore using $Results/agescore.tex, ci(3) b(3)  replace label 
	
			  
			  
		collapse (sum) Obs = ones , by(agebins)
		
		label variable Obs "Observations (000s)"
	
	replace Obs = Obs/1000
		
	
	coefplot 	(agescore, keep(1.agebins) ciopts(recast(rcap) lcolor(red)) msymbol(D) mcolor(red) msize(large) ) ///
	(agescore, keep(2.agebins) ciopts(recast(rcap) lcolor(orange))  msymbol(D) mcolor(orange) msize(large)) ///
	(agescore, keep(3.agebins) ciopts(recast(rcap) lcolor(gold)) msymbol(D) mcolor(gold) msize(large)) ///
	(agescore, keep(4.agebins) ciopts(recast(rcap) lcolor(green)) msymbol(D) mcolor(green) msize(large)) ///
	(agescore, keep(5.agebins) ciopts(recast(rcap) lcolor(blue)) msymbol(D) mcolor(blue) msize(large)) ///
	( agescore, keep(6.agebins) ciopts(recast(rcap) lcolor("111 000 255")) mcolor("111 000 255") msymbol(D) msize(large) ) ///	
	, vertical levels(95) yline(0)  nooffsets  legend(off) xlabel(, angle(60))  msymbol(D) omitted baselevels xscale(range(1 6)) ///
	yaxis(2) ytitle(Difference in score relative to < 30, axis(2)) coeflabels( ///
                1.agebins = "<30" ///
			  2.agebins = "30" ///
			   3.agebins = "40" ///
			  4.agebins = "50" ///
			  5.agebins = "60" ///
			  6.agebins = "{&ge}70" ///
			  )
	
	
	
	
			
				

		
		addplot: bar Obs agebins if agebins ==1, yaxis(1) bcolor(red%20) xscale(range(1 7))  xlabel(, valuelabel) xtick(7)
		addplot: bar Obs agebins if agebins ==2, yaxis(1) bcolor(orange%20) xscale(range(1 7)) xlabel(, valuelabel) xtick(7)
		addplot: bar Obs agebins if agebins ==3, yaxis(1) bcolor(gold%20) xscale(range(1 7)) xlabel(, valuelabel) xtick(7)
		addplot: bar Obs agebins if agebins ==4, yaxis(1) bcolor(green%20) xscale(range(1 7)) xlabel(, valuelabel) xtick(7)
		addplot: bar Obs agebins if agebins ==5, yaxis(1) bcolor(blue%20) xscale(range(1 7)) xlabel(, valuelabel) xtick(7)
		addplot: bar Obs agebins if agebins ==6, yaxis(1) bcolor("111 000 255 %20") xscale(range(1 7)) xlabel(, valuelabel) xtick(7) ytitle(Observations (000s), axis(1))
		
		
		
		
			graph export $Results/agescore_obs.png, replace

			
			
			
			
		
		
			
		estpost tabstat Obs , by(agebins) 	statistics(mean ) columns(statistics) listwise

 
		
		esttab . using $Results/ageobs.tex, main(mean) nostar unstack noobs nonote label mtitle("1000 Observations by Age") replace
		
